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The relaxation dynamics of correlated electron transport (ET) along molecular chains is studied 
based on a substantially improved numerically exact path integral Monte Carlo (PIMC) approach. 
As archetypical model we consider a Hubbard chain containing two interacting electrons coupled to 
a bosonic bath. For this generalization of the ubiquitous spin-boson model, the intricate interde- 
pendence of correlations and dissipation leads to non-Boltzmann thermal equilibrium distributions 
for many-body states. By mapping the multi-particle dynamics onto an isomorphic single particle 
motion this phenomenon is shown to be sensitive to the particle statistics and due to its robustness 
allows for new control schemes in designed quantum aggregates. 
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Recent progress in engineering and manufacturing of 
nanoscopic electronic devices enabled to tackle the long 
standing problem of contacting single molecules P, 0, U 
This remarkable achievement not only opens a way for 
building the ultimate transistor on the molecular basis, 
but is also expected to become an important tool for 
investigating electron transfer (ET), which is one of the 
most fundamental processes in chemistry and biology , 
the best known example being photosynthesis 0. Like- 
wise, in semiconductors charge transfer plays a key role 
and most recently studied arrays of few electron quan- 
tum dots may be seen as artificial molecular structures 
which, in contrast to their native counterparts, allow for 
a completely controllable ET . 

The most successful paradigm for the description of 
ET is the donor-acceptor model. In its simplest form it 
is represented by a two-site chain where a single electron 
is transferred between the sites via tunneling hybridiza- 
tion. However, irreversible transfer can only take place 
in presence of an environment consisting e.g. of residual 
vibronic degrees of freedom or electromagnetic fluctua- 
tions. In the past, the Caldeira-Leggett model has 
been shown to capture the relevant features of dissipative 
ET, where the electronic system is coupled to a bath of 
harmonic oscillators . In case of only two sites and one 
single charge this is the famous spin-boson model, which, 
together with its generalizations to multi-site geometry, 
is a paradigmatic model on its own comprising e.g. quan- 
tum Brownian motion, Kondo physics, Luttinger liquids 
and atomic quantum dots Despite its simplicity 
the associated phenomenology is extremely rich, showing 
different sorts of quantum phase transitions [LU^ . 

With the advent of fully tunable quantum structures a 
fundamental constraint has to be relaxed though, namely, 
the restriction to single particle dynamics. Multi-particle 
correlation effects may undoubtedly have a potential to 
considerably alter the simple ET physics and open the 
door for new control mechanisms ^M- Of course, 
correlated ET in a dissipative reservoir is a formidable 
task, mainly due to the fact, that a quantum reservoir 
leads to strong retardation effects. So far virtually all 
studies made use of perturbative approaches focusing 



on either completely coherent transfer with no dissipa- 
tion |l3l Il4j , or completely sequential transfer with very 
strong Coulomb repulsion |l5l |. which effectively reduces 
to a single particle ET problem. The inclusion of both 
correlations and dissipation has been attacked only re- 
cently in |l6| |. where equilibrium properties of a two site 
spinful system filled with two electrons has been analyzed 
with the help of numerical renormalization group (RG) 
technique. However, being intrinsically an equilibrium 
technique, this method cannot reveal any information 
about dynamical properties of these systems. In addi- 
tion, finite temperature effects are usually quite difficult 
to access in the framework of the RG. In this Letter we 
close this gap and, based on a numerically exact path 
integral Monte Carlo (PIMC) scheme, we present results 
for the real-time dynamics of two correlated fermions in 
a dissipative environment. This model can be seen as the 
simplest generalization of the single particle spin-boson 
model to include many-body effects. It may thus reveal 
archetypical properties of correlated dissipative quantum 
dynamics. 

The system is modeled by an open Hubbard chain with 
N sites of spacing a 



N 



^ A, (4 + /i.e.) , (1) 
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where di^ are annihilation operators for electrons with 
spin a on the site i. Ei are the bare energies of the 
levels, Ui the corresponding interaction strengths and 
Ai are the tunneling matrix elements. The interac- 
tion with the bosonic bath given by the Hamiltonian 
Hb — J2a (-fa/^^^Q + maUj'^Xa/2) IS accomplished via 
a standard dipole coupling 



Hi = -aV CaXa + a^V^ 



(2) 
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where 



N 



V=J2 [*-(^ + l)/2]rfU. 



(3) 



is the polarization operator of the Hubbard chain. 

We are interested in the reduced dynamics of the 
fermionic system determined by the density opera- 
tor p{t) = TrB{cxp{-iHt/h)W{0) exp{iHt/h)}, where 
W{0) specifies the initial state of the compound. Here 
the trace over the bosonic bath degrees of freedom can be 
done exactly within the path integral approach by rep- 
resenting the density operator in the site representation 
of the chain in terms of a pseudo-spin- (iV — l)/2 opera- 
tor S"' . Accordingly, the eigenvalues of S'f represent the 
position of the respective electron along the chain, i.e. 
S^\s'') = s'^ls'^) (s'^ = -S,...,S), and Is"") denote a lo- 
calized single-electron state. In the rest of the Letter we 
focus on a chain filled with two electrons of opposite spins 
with a = 1,4. Then, the corresponding many-body basis 
{|s^,s^)} follows from properly anti-symmetrized single 
particle states. The relevant dynamical observables are 
the populations P\{t) = Trs{^p{t)} where A denotes a 
projection operator either onto a certain site and spin or 
onto a certain many-body state. Initially, the system is 
prepared in a specific many-body state with the bosonic 
reservoir equilibrated to that state. 

This way, the reduced density operator is expressed as 
a double path integral along a Keldysh contour with for- 
ward s'^ and backward paths. The impact of the dis- 
sipative environment appears as an influence functional 
introducing arbitrarily long-ranged interactions in time 
between the paths. It is convenient to switch to the com- 
binations jy'^ = 3°' + s"' and — s"' ~ s"' so that one 
arrives at the exact expression 

PaW = /l?r;/l?^^A(77,0^[^,| exp(-$[,7,|) , (4) 



where ff = {rf,?]^) and ^ = (C^,^'^). Here, A is the bare 
action factor in absence of a reservoir and the influence 
functional reads 
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with e = (1, 1). The kernel L{t) = L'{t) + iL"{t) is re- 
lated to the force-force auto-correlation function of the 
bath and is completely determined by the spectral den- 
sity of its modes. Further, p. = lim;^^_>o S/3i(0). The 
immediate advantage of the representation Q is that 
formally it looks like the dynamics of a single particle in 
a two-dimensional plane. We shall discuss this point in 
more detail below. Note further, that even in absence of 



a direct Coulomb interaction the two charges are corre- 
lated due to the coupling to the heat bath. 

An analytical treatment of the expression is in gen- 
eral not feasible, mainly due to the retardation in the 
influence functional which grows with decreasing tem- 
perature. In this situation PIMC methods have been 
shown to be very powerful numerically exact means to 
explore the non-perturbative range. Our starting point 
is a PIMC scheme ^3 which exploits the linear depen- 
dence of ^'[?f, on the quasi-classical paths ^(i), so that 
in (0J) the corresponding summations can be expressed 
as a scries of simple matrix multiplications and there- 
fore be carried out explicitly. These matrices include 
the bare short time propagators of the system and phase 
factors stemming from the rf-dcpendent part of <&. Con- 
sequently, the number of MC variables is reduced by a 
factor of two, which leads to a significant soothing of the 
dynamical sign problem |l8j . Previous simulations were 
restricted to single particle transport on up to three sites 
and for shorter times only, because the standard MC- 
weight requires to perform the full scries of matrix multi- 
plications for each MC-stcp due to the retardation of the 
rf-dependent phase factor. A breakthrough is gained by 
neglecting the long-time retardations during the propaga- 
tion of the MC trajectory. Since they are fully accounted 
for during the final accumulation process, the numerical 
exactness of the MC scheme is not impaired. This way, 
one achieves a strong decoupling of quasi-classical (ff) and 
quantum (^) coordinates, which allows to calculate and 
store products of short time propagators independent of 
the MC-sampling. Eventually a speed-up with respect 
to the original method ^3 by a factor of about 100 is 
gained. 

The simplest non-trivial system which is able to accom- 
modate two electrons and thus the simplest many-body 
generalization of the spin-boson model is a two-site chain. 
The proper representation of a corresponding many-body 
basis is given by 



s = {UT,o), U,T), IT,i), |o,it>}, 



(5) 



where the left (right) argument of the ket represents the 
spin occupation of site — l/2(-|-l/2). For these localized 
spin states the relevant observables are the site popu- 
lations and the populations occupying one of the 
states e B. In particular, the latter ones give de- 
tailed insight into the intimate relation between correla- 
tions and dissipation during the relaxation from an initial 
non-equilibrium state. For very long times, the expec- 
tation is that the system reaches a thermal equilibrium, 
where all four many-body states are populated according 
to a Boltzmann distribution. Thus, in case of vanishing 
Coulomb interaction Ui = and energetically completely 
degenerated sites Ei = 0, one anticipates, no matter how 
the system is prepared at t = 0, a final probability of 
1/4 for each state in Q, whereas both sites of the chain 
have to be half-filled. While the latter expectation turns 
out to be true, remarkably, the former one cannot be 
confirmed by our simulations (see Fig. Not only do 
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FIG. 1: Relaxation dynamics of two correlated electrons in a 
degenerate two-site systems according to Fig. Shown are 
the populations for the states | J,|, 0) (circles), | J,, t) (triangles 
up), I Ti i) (triangles down), and |0,J,t) (diamonds) for a bath 
with ohmic spectral density and exponential cutoff at uJc = 
5 A; the damping strength is a = 1/4 and temperature /1/3A = 
0.1. Dashed lines represent theoretical equilibrium values. 
Inset: Dynamics for a system according to Fig. Error 
bars (vertical lines) are typically of the size of the symbols. 



the populations Pip{t) saturate at different values, but 
these also depend on the initial preparation. This at 
first startling dynamics leading to a non-Boltzmann dis- 
tributed steady state in the many-body basis |SJ), can 
be understood in terms of decoherence free states (DFS) 
[T^ . These states arc simultaneous eigenstates of Hg 
and of the polarization operator V. Accordingly, after 
being prepared in such a state the system stays there 
forever and no relaxation into the equilibrium can occur. 
While DFS have been studied for very weakly damped 
systems in the context of quantum information process- 
ing, where the generic basis is the eigenstate basis of Hs, 
in strongly condensed phase systems the site represen- 
tation is the proper one. Hence, for correlated ET the 
impact of DFS onto transport properties is important and 
leads to completely new and unexpected phenomena. In 
the case considered here, there is only one such state, 
namely 

l<^i) = (U,T>-|Ta»/x/2, 

which corresponds to zero dipolc moment and apparently 
is a delocalized spin state. An alternative basis set of the 
electronic Hilbert space is thus given by 
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l0i),^(U,T) 



IT,i)),UT,o), |o,iT) 



consisting of both delocalized {\4>i), 102)) and localized 
(103), 104)) spin states. Obviously, |0i) does not partic- 
ipate in the relaxation process, such that equilibration 
is restricted to the sub-space spanned by |0i), I = 2,3, 4, 
leading for an energetically degenerate system to 

P^,(t^oo) = [1-P^,(t = 0)]/3 , 1 = 2,3,4 (6) 
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FIG. 2: Two-site Hubbard chain with standard (a) and spin 
flip/two particle tunneling (b) couplings. 



Now, due to Piifi 
+ P^, . and P? 



for t - 



= P^„ Pu + P^.l = 

"I" -f</)i J anu r-^--^ = P|°^ for t oo due to symmetry, 
one easily understands from ^ both the non-Boltzmann 
distribution of the localized states as well as their depen- 
dence on the initial preparation. Imagine that initially 
the electronic sub-system is prepared in one of the states 
in jSJ, which experimentally is the typical situation. 
Then, for the non-equilibrium dynamics two different sce- 
narios are possible: If both electrons initially occupy the 
same site, i.e. -PxT,o(0) = 1 or Po,iT(0) = 1, one always 
has (0) = so that the final equilibrium distributions 



read Pi 



P 



o,iT 



1/3 and P 



poo 



1/6. If, 



how- 
ever, both sites are initially occupied, i.e. P|_|(0) = 1 or 
P|,x(0) = 1, so that P01 (0) = 1/2, the final distributions 
become P|^o = P;^^^ = 1/6 and P~ = Pf°^ = 1/3. 
Therefore, the system, after having been prepared in 
one of the states always remembers, whether this 
state had a vanishing or non- vanishing dipolc moment. 
This phenomenon crucially depends on the symmetry of 
the operator mediating the interaction with the bosonic 
baths. However, since the dipole coupling specified in 
(0) is generic for a variety of native and artificial molec- 
ular structures, the above dynamics can be expected to 
be typical. In most native molecules ET is triggered by 
photo-excitation from an electronic ground state so 
that correlated ET is difficult to observe. In contrast, 
e.g. in synthetic aggre gate s, which are very weakly cou- 
pled to external leads [Tj, [i3 , or in arrays of quantum 
dots , the excess charge can be better controlled. Then, 
starting e.g. from a totally unpolarized ensemble, after 
equilibration only 1 /3 of its constituents would be polar- 
ized and 2/3 of them mipolarizcd. This is an example of 
a non-Boltzmann type of thermal distribution, which for 
the generic form of system-bath coupling is a signature 
of: (i) the topology of the tunneling couplings in Hs [cf. 
(^3)] and (ii) Fcrmi-Dirac statistics. Let us analyze these 
two points in more detail. 

To address (i), one first realizes that the existence of a 
common eigenstate |(/)i) of Hq and V is completely in- 
dependent of the electronic energies Ei, the Coulomb 
interaction strengths Ui, temperature and the spectral 
density of the bath modes. Accordingly, the described 
phenomena can also be observed for non-degenerate sys- 
tems (El 7^ i?2 7^ and/or Ui =^ U2 0)- In this case. 
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FIG. 3: Same as in Fig. [H but for hf3A = 0.5. Inset: Dy- 
namics for a system with ft/3 A = 0.1, Coulomb interaction 
Ui = U2 = 4ftA, and Ei = E2 = 0, i.e. (iT,0|i^s| iT,0) = 
(0, iT l^^slO, it) = 4ftA, (i, T \Hs\ I, T) = (T, i \Hs\ T, i> = 0. 
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FIG. 4: Same as in Fig. but for two identical bosons on 
a chain with A'^ = 3. Shown are the populations for states 
|n_i,no,ni) with occupation m on site i: |2,0, 0) (circles), 
|1,1,0) (squares), |1,0, 1) (triangles up), jO, 2,0) (triangles 
down), |0, 1, 1) (diamonds), and |0,0, 2) (triangles right). In- 
set: Dynamics for N = 2 with states |?i-i/2, 71-1/2); |2,0) (cir- 
cles), 1 1, 1) (squares), jO, 2) (triangles). 

one even arrives at the counter-intuitive situation that 
in thermal equilibrium the occupations of energetically 
higher lying states exceed those of energetically lower ly- 
ing ones (cf. inset in Fig. (S)). Moreover, effects tied to 
the DFS also survive the coherent-incoherent transition, 
its typical signature being the onset of oscillations in the 



population dynamics (cf. Fig. Thus, they prevail in 
two distinct dynamical regimes. However, DFS crucially 
depend on the topology of the tunneling couplings in Hs- 
They always exist for isotropic couplings, but can be eas- 
ily destroyed e.g. by appropriately adding or removing 
couplings between the localized states (0. For exam- 
ple, interchanging the coupling according to Fig. so 
that spin flips and two-particle tunneling arc allowed, the 
Boltzmann equilibrium is restored (sec inset in Fig. 

To discuss (ii), one easily deduces a mapping of the dis- 
sipativc dynamics of two correlated fermions on two sites 
onto a single particle motion along the edges of a two- 
dimensional square, spanned by the many body states 
(O , see Fig. This mapping has two immediate conse- 
quences that will be studied in details elsewhere: First, 
it allows to adapt perturbative approaches like the NIBA 
and the NICA |3| developed for ordinary spin-boson mod- 
els to capture also correlated many-body dynamics; on 
the other hand, the dynamics along any A^-site chain with 
an arbitrary number of electrons can be mapped onto 
a dissipative single particle ET on a higher-dimensional 
square For indistinguishable bosons such a map- 

ping is also found leading though to a triangular lattice. 
Then one derives from the geometry of the mapping that 
for a Hubbard chain only one DFS exists whereas for 
the similar bosonic model such states are only possible 
for odd N. This is illustrated in Fig. 01 for chains with 
A = 3 and N = 2. In the former case one has final pop- 
ulations of 1/5 vs. 1/10, while in the latter one they all 
saturate at 1/3. Note that the A^ = 3 fcrmionic system 
exhibits smaller equilibrium values of 1/8 vs. 1/12 due 
to the larger number of accessible many-body states. Of 
course, DFS never exist in one-dimensional single parti- 
cle systems. 

To summarize, we have developed a PIMC technique 
which opens the door for a reliable modeling of the 
nonequilibrium dynamics of correlated ET even on longer 
time scales. As a remarkable many-body effect, related 
to the existence of a DFS, a non-Boltzmann steady state 
was predicted. A mapping onto an isomorphic single par- 
ticle model reveals its sensitivity to the particle statistics. 
Experimentally, this phenomenon may lead to new con- 
trol schemes for fermionic or bosonic transport in tight- 
binding systems like synthesized molecular chains, quan- 
tum dot arrays, or optical lattices. 
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